Thermal recoil force, telemetry, and the Pioneer anomaly 
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Precision navigation of spacecraft requires accurate knowledge of small forces, including the recoil 
force due to anisotropics of thermal radiation emitted by spacecraft systems. We develop a formalism 
to derive the thermal recoil force from the basic principles of radiative heat exchange and energy- 
momentum conservation. The thermal power emitted by the spacecraft can be computed from 
engineering data obtained from flight telemetry, which yields a practical approach to incorporate 
the thermal recoil force into precision spacecraft navigation. Alternatively, orbit determination can 
be used to estimate the contribution of the thermal recoil force. We apply this approach to the 
Pioneer anomaly using a simulated Pioneer 10 Doppler data set. 

PACS numbers: 95.10.Eg,95.30.Sf,95.55.-n,95.55.Pe 



I. INTRODUCTION 

Precision navigation of spacecraft requires a detailed 
knowledge of small forces affecting its trajectory. This 
includes the recoil force associated with thermal radia- 
tion. It recently became clear that the small force due 
to the radiation of waste heat from the spacecraft itself 
cannot be ignored [1.]. This is especially true in the case 
of the Pioneer 10 and 11 spacecraft. A small, anomalous 
Doppler residual is present in these spacecrafts' radio sig- 
nal, and it is believed to be caused by an acceleration of 
unknown origin [3, [3, 3 ■ Investigations of this anomaly 
indicated that its magnitude is comparable to that of the 
acceleration due to a thermal recoil force 0, H, Q . There- 
fore, the possibility that the anomalous acceleration is 
caused by an anisotropic emission of the heat generated 
on board of the spacecraft must be investigated 

Surprisingly, there is no established formalism for in- 
corporating the thermal recoil force into precision space- 
craft navigation^. In fact, most of the literature is deal- 
ing with energy distribution processes, but not momen- 
tum transfer. As a result, some authors were forced to 
develop ad hoc formalisms in order to estimate the ef- 
fects of the thermal recoil force on spacecraft trajectories 
@, Sin [11; while others (e.g., [|) only roughly esti- 
mate the magnitude of this force, but do not precisely 
determine its direction and temporal dependence, nor do 
they provide formal error estimates. However, this is 
exactly what we need to investigate the nature of the 
Pioneer anomaly. 

For navigational purposes, one is interested in the mag- 



*URL: 'http : //www . vttoth . com/J_ 



nitude, direction, and temporal evolution of forces act- 
ing on a spacecraft [H. For the recoil force due to 
internally generated heat these can be estimated using 
available knowledge of the spacecraft's geometry, ther- 
mal properties, and state. The geometry and thermal 
properties of the spacecraft (excluding effects of aging) 
can be established from design documentation, prelaunch 
test results, and calibration experiments performed after 
launch. The spacecraft's state, in turn, can be obtained 
for any given moment of time from its telemetry, which 
contains unique data on the power consumption, thermal 
status, and physical configuration of the craft [l|]. 

The purpose of this paper is to develop the methods, 
tools, and procedures that are needed for the evaluation 
of the thermal recoil force as a likely cause of the anoma- 
lous acceleration of the Pioneer 10 and 11 spacecraft. 
This investigation of the on-board small forces is per- 
formed in conjunction with the analysis of the Doppler 
data, allowing us to show, for the first time, how preci- 
sion orbit determination can also be used to estimate a 
spacecraft's thermal properties. 

This paper is organized as follows. Our analysis be- 
gins with considering heat conduction and radiation in 
Secini We calculate radiation pressure in Sec. lIIII Next, 
in Sec. IIVI we explore in detail the resulting recoil force 
and its relationship with the amount of heat generated. 
In Sec. |Vl we analyze the accuracy with which the recoil 
force can be determined. Next, in Sec. IVIl we connect 
the computation of the recoil force with orbital analy- 
sis. In Sec. IVIIl we apply what we learned to simulated 
Pioneer 10 orbital data. We present our conclusions in 
SecEml 
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II. HEAT CONDUCTION AND RADIATION 

Heat flows from internal heat sources to external ra- 
diating surfaces via three mechanisms: convection, con- 
duction, and radiation. For the purposes of this paper, 
we ignore convection within the spacecraft^. 

Heat conduction is described by Fourier's law 



-fcVT, 



(1) 



where q is the heat flux (measured in units of power over 
area), T is the temperature, and k is the heat conduction 
coefficient. In the general case, fc is a tensorial quantity, 
but for homogeneous and isotropic materials, k reduces 
to a scalar coefficient^. In general, q, fc, and T are all 
functions of the coordinates x and time t. 

Heat flux also obeys the energy conservation equation 



V-q = 6-Cft.p— , 



(2) 



where b is the volumetric heat release (measured in units 
of power density), Ch is the material's specific heat, and 
p is its density. 

Equation ^ has general applicability. In many cases 
of practical interest the internal heat sources are com- 
pact and pointlike, and their thermal power is known: 
for instance, instrumentation boxes within a spacecraft 
compartment. Denoting the thermal power of n thermal 
sources as Bi(t) {i — l...n), we can write 



fe(x,t) 



i=l 



(3) 



where is the location of the ith heat source and 5 is 
Dirac's delta function. 

A system is in a steady state if its properties do not 
change with time. In particular, this means dT / dt — 0, 
which leaves us with 



V • q = 6. 



(4) 



The heat conducted to a surface element must equal 
the heat radiated by that surface element. Therefore, at 
the surface. 



<? = q • a. 



(5) 



where a is the unit normal of the radiating surface ele- 
ment, and q is the surface element's radiant intensity. 



We note that convection may be significant in the case of space- 
craft that utihze a liquid or gas coohng system or in which sub- 
stantial quantities of fuel can flow from one part of the spacecraft 
to another. 

A notable case of anisotropic conductivity where k is tensorial is 
that of multilayer insulation. 



The radiant intensity, or energy flux, of a radiat- 
ing surface is related to its temperature by the Stefan- 
Boltzmann law: 



g(x,<) = ae(x,t,r)r4(x,i), 



(6) 



where cr ~ 5.67 x lO"* Wm^^K^^ jg ^Yie Stefan- 
Boltzmann constant, while the dimensionless coefficient 
< e < 1 is a physical characteristic of the emitting sur- 
face. This coefficient can vary not only as a function of 
location and time, but also as a function of temperature. 

We describe radiation by a four-dimensional stress- 
energy-momentum tensor that takes the form 



c '^u p 

p p 



(7) 



where u is the energy density of the radiation field, p is its 
momentum density, and P is the radiation pressure ten- 
sor. {T^^'^ is, in fact, the stress-energy-momentum tensor 
of the electromagnetic field in a vacuum.) 

The stress-energy-momentum tensor obeys the conser- 
vation equation 



(8) 



where denotes covariant differentiation with respect 
to the coordinate . In three dimensions, this yields the 
following conservation equations: 



c 



Vp 

V -P 



0, 
0, 



(9) 
(10) 



where the dot denotes differentiation with respect to 
time, i.e., x = dx/dt. 

The energy E' in a given volume V is 



E = 



idV 



by definition. The power, denoted by Q, is 



Q = ^ ^ [udV = c 
dt J 

or, after applying Gauss's theorem. 



Q 



pdA 



V • pdV, 



p ■ adA, 



(11) 



(12) 



(13) 



where A represents a surface enclosing the volume V and 
a is the unit normal of surface element dA with area dA, 
such that dA = adA. 

Comparing with ([5|) and noting that Q = J qdA, we 
obtain the relationship between radiant intensity and mo- 
mentum density at the radiating surface: 



q(x,i) = c2p(x,i) • a. 



(14) 



We describe the radiative flow of energy E using the 
intensity^ /, which is the flow of energy across surface 



* Some textbooks call the quantity / the radiance, and its integral 
over a finite surface the (radiant) intensity. 
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element dA, in a time interval dt, in the solid angle duj 
around direction n: 



E = 



/(x, t, n)n • dAdujdt, 



or, after differentiating with respect to t 
dE 



(15) 



Q = :i- = /(x, i, n)n • dAdw. (16) 

Comparing with , we get 

p = ^ y" /(x,i,n)ndw. (17) 

From p4)) . then, we obtain 

q{x,t) = J I{x,t,n)n- aduj. (18) 

A surface is defined as a diffuse or Lambertian emitter 
if the intensity / does not depend on the direction of 
radiation emanating from a surface element. In this case, 
we can take / outside the integral sign and write 



q{x,t) = I{x,t) J n-aduj. 



(19) 



The integral on the right-hand side should be evalu- 
ated over a hemispherical surface of unit radius centered 
around the surface element dA. Parameterizing the inte- 
gration surface using spherical coordinates (r, (f), 9) (with 
= at the north pole), we note that du = sm (j)d(j)d9 
and n • a cos 0, and the integral reads 

n • adhj = / cos (j) am(j)d(j)d6 — tt, (20) 
Jo Jo 

thus 

q{x,t)^TTl{x,t). (21) 

Radiation from a surface element dA = dAi that is 
intercepted by a second surface element dA2 at distance 
r can be calculated by using, as the solid angle, dw = 
r-^n ■ dA2 in ^6]) : 

/(x,t,n) 



11^2 



-n • dAin • dA2 



/(x, t, n) cos 9i cos 92 



dAidA2, (22) 



where 9i and 6*2, defined by cos 6*1 = n ■ dAi/dAi and 
cos 02 — n • dA2/dA2, are the angles between the di- 
rection of heat radiation and the normals of the surface 
elements dAi and dA2, respectively. 

If both surfaces are Lambertian emitters, and we take 
into account heat flowing in both directions, the heat Q12 
exchanged between the two surfaces is 



gr2 = / / fa-'^^)^"^^^ieos02 ^^^^^^ 



TTr 

(eiTf - e2T^)cos9i COS6I2 



These results describe the radiative exchange of energy. 
Next, we turn our attention to momentum exchange. 



III. RADIATION PRESSURE AND THE 
RECOIL FORCE 

The pressure tensor of radiation is written as 

P(x, ^(X' n)nndu;, (24) 

where uv is the dyadic product of two vectors u and 
V. (This form of the pressure tensor coincides with the 
Maxwell stress tensor for plane or spherical electromag- 
netic waves.) For a Lambertian emitter, once again we 
can take / outside the integral sign, yielding 



''(x,t) = -I{x,t) j nnduj. 



(25) 



The recoil force acting on the emitter will be of the 
same magnitude, but opposite in sign to the change in 
momentum in a given volume. It can be written as 



F{t)^-J pdV. 



(26) 



After using ^TU^ and then applying Gauss's theorem, we 
get 

F(t) = - ■ = - y P(x, t) ■ dA. (27) 

For a Lambertian emitter, this force then becomes [gl-ITol. 



F{t) = J I{x,t) J nn- aduj 



dA. 



(28) 



The inner integral can be evaluated by making use of the 
identity (ab) • c = a(b • c): 



J (nn) • aduj = J {a - n)nduj. 



(29) 



As with (|19p. we integrate over a hemispherical surface 
of unit radius centered around the surface element dA. 
We set up a spherical coordinate system (r, cj), 9) such that 
= corresponds to the north pole (that is, the direction 
of the surface normal a), and we set up two additional 
basis vectors b(0 — 7r/2, 9 ~ Q) and c(0 — 7r/2, 9 = tt/2) 
such that n can be expressed as 

n = (n • a)a + (n • b)b + (n • c)c 

= cos (pa + sin (j> cos 6h + sin (j> sin 9c. (30) 

To integrate , we use a • n = cos <j) to obtain 



dAidA2. (23) J cos0[cos(?!)a + sin(/>cos6'b + sint/isin^cjdti 



2" 2 2^ 

/ cos (j) sin (j)ad(^d9 — —a. 
Jo 3 



(31) 
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Therefore, 



F{t) = -^^ I g(x,t)dA, 



(32) 



indicating that, in the Lambertian case, the radiation 
pressure is isotropic, and the radiation pressure tensor 
reduces to a scalar quantity: 



3 c 



(33) 



where I is the identity tensor. 

As described previously, q can be obtained by solving 
the heat conduction equations (H]) and along with 
the radiative heat transfer equation (|23p and boundary 
conditions. 

The surface density of the recoil force that corresponds 
to p3p . acting on surface element dA with unit normal 
a, is 



2 1 

f(x,<) = ---q{x,t)a.. 

3 c 



(34) 



From this, the recoil force can be computed by integra- 
tion^ : 

F(t) = J f(x,t)dA = -^i J g(x,t)dA. (35) 

This result establishes the relationship we sought be- 
tween radiated heat and the associated recoil force. 



IV. THE RECOIL FORCE AND ITS SOURCES 

Heat conduction q inside a heat emitting object and 
the radiant intensity q at its exterior surfaces can be ob- 
tained by solving, in the steady state, Eqs. ([4]) and ^ 
for the unknown functions q and q using the known func- 
tions a (representing the emitter's geometry) and b (the 
volumetric heat release inside the emitter), along with 
appropriate boundary conditions (e.g., sky temperature). 
Given constant boundary conditions and an unchanging 
geometry, the solution for (/(x, t) can be obtained in terms 
of 6(x, t). The recoil force F{t), which is a functional of 
g(x,t) as per ([35)1 . can therefore be expressed as a func- 
tional of 6(x, t): 



F(0=F[6(x,i)]. 



(36) 



^ Knowing the recoil force surface density also allows us to compute 
the torque acting on the emitter. The torque surface density is 
X = (x — xq) X f, where xq is the location of the the emitter's 
center-of-gravity. The total torque T, then, can be written as 



[-(t) = J (x-xo) X {{x,t)dA = y 9(x,t)(x - xo) 



X dA. 



In many cases, b can be represented by discrete, compact 
heat sources in accordance with ([3|). As we noted, this 
is the case, in particular for spacecraft containing instru- 
mentation boxes within its compartment, with teleme- 
tered power readings available for each. In this case, we 
can write F as a functional of the n functions Bi{t): 



(37) 



The magnitude and direction of the recoil force are 
both functions of time. However, if the emitting object 
is rotating, and its rate of rotation is sufficiently high 
compared to the rate of change of the recoil force, the 
time average vector components of the recoil force that 
lie in the plane of rotation will be negligible. The residual 
recoil force will always be perpendicular to the plane of 
rotation, i.e., parallel to the rotating object's spin axis. 

To see this, we write the recoil force in the form 



where 



F(t) = F||(t)+M(L^t)-Fi(0, 



F||(t) = (F(t).s)s 



(38) 



(39) 



is the component of F{t) parallel with the spin axis rep- 
resented by the unit vector s (which, we assume, remains 
constant in time), and 



F^{t)^R-\Luto){F{t)~F\^{t)) 



(40) 



is a perpendicular component of F(t) in a corotating ref- 
erence frame, while K((/)) is a tensor representing a rota- 
tion by the angle (j). 

Ignoring forces other than the recoil force, the posi- 
tion of the rotating object as a function of time can be 
calculated as 

x(t) = x(to) + Mto)it -to)+ J J ^Fit)dh, (41) 

where x(io) is the position, x(to) the velocity of the ob- 
ject at some time t — tQ. We assume that the object's 
mass, m, remains constant in time. We denote the dis- 
placement of the object as a result of the recoil force as 
Ax = x(t) — Xq — Voi. Therefore, 



Ax=— J J {F\\{t) + M.{ujt) ■Fj_{t))d^t. 



(42) 



The displacement due to F|| is entirely along the spin 
axis. To calculate the displacement due to the perpendic- 
ular component, we assume that it is well approximated 
by a linear function of time: 



F±{t) 



Fit, 



(43) 



where F^ is constant. Thereafter, noting that 
J M.{LLjt)dt — uj~^'R{ujt — 7r/2), we calculate the dou- 
ble integral for n full revolutions over a time period 
At = 27m/w and obtain 
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Ax^ = -^r—R(LutQ - tt) • F_LAt, (44) 

describing an arithmetic spiral. In three dimensions, 
then, the motion of the object is described by a heUx of 
widening radius around the object's original axis of ro- 
tation. The growth of the width of the helix is governed 
by F^. If Fj^ remains nearly constant in time, F^ ~ 
and the object's spin axis returns to its original position 
after every full revolution with no cumulative displace- 
ment in the perpendicular direction. Therefore, ignoring 
a constant F^ introduces only a small periodic error and 
no cumulative error in the calculations, and ignoring Fj^ 
in the case when F_l is nonvanishing but small and ap- 
proximately constant only introduces a small cumulative 
error. This means that in the case of a spinning object, 
the recoil force (|37p can be expressed as 

F(i) -F||[i?i(t),B2(0,-,Sn(i)]s- (45) 

Without loss of generality, F|| can be expanded in the 
form of a Taylor series in the Bi. For the purposes of the 
present paper, it is sufficient to keep only linear terms. 
With this in mind, and noting that -F||(0, ...,0) = 0, we 
can write 

1 " 

F\\it) = -J2C^B,it), (46) 
where the dimensionless coefficients are given by 

The factors are determined by the geometry and opti- 
cal properties of the emitter, and are expected to remain 
constant so long as the emitter's geometry and optical 
properties do not change. 

For any given heat source, the principle of conserva- 
tion of energy dictates that — 1 < < 1. The coeffi- 
cient is zero if heat from that particular source is emitted 
isotropically, resulting in no net recoil force. Therefore, 
these factors determine the efficiency with which each 
heat source contributes to the object's acceleration. 

V. APPLICATION AND ERROR ANALYSIS 

The formalism that we obtained can be employed in 
a direct calculation of the recoil force using conventional 
numerical methods of heat transfer. Together, Eqs. ([!]) 
and ([2]) are two first-order differential equations in two 
unknown functions q and T that describe heat conduc- 
tion inside materials, while represents additional 
constraint equations describing radiatively coupled sur- 
face elements. 

A specific solution can be obtained if the material prop- 
erties (represented by k and e, as well as Ch and p) and 



heat sources (represented by h) are known, and appropri- 
ate boundary conditions are given. 

One such boundary condition is the steady-state con- 
dition dT/dt = 0, in which case we replace © with 
For a heat emitting object situated in empty space, an- 
other boundary condition can be specified in the form 
of the sky temperature (i.e., the microwave background 
radiation temperature) to which the object's exterior sur- 
faces are radiatively coupled. A practical difficulty arises 
if the object has facing surfaces (i.e., surfaces that are 
radiatively coupled to one another) but these situations 
are dealt with easily using standard finite element or ray- 
tracing numerical codes. 

It should be noted that in this solution is 

fully specified when the volumetric heat release 5(x, t) 
is known, even as no temperature values inside the emit- 
ting object are given. When spacecraft telemetry pro- 
vides both power and temperature measurements, these 
data together represent a redundant data set that can 
be used to verify and validate thermal models. Labori- 
ous but, in principle, straightforward application of these 
equations can lead to a temperature map of the exterior 
surfaces of an emitter. When this temperature map is 
known, the recoil force can be computed directly using 
Eq. dSSl). 

An alternative to evaluating the vector-valued integral 
psp along the nontrivial exterior geometry of the emit- 
ter is the use of a control volume technique iTj. The 
anisotropy of thermal emissions can be determined by 
surrounding the object with an infinite control volume, 
which is approximated by a sufficiently large fictitious 
spherical surface centered around the emitter that is used 
to intercept all radiated heat coming from the emitter 
(see discussion of a similar approach involving pixel ar- 
rays in [l3|). The recoil force is computed by evaluating 
the integral ([35|) along the surface of this sphere. These 
calculations can be carried out using standard thermal 
modeling tools^, that have been used successfully for 
the design and operations of many missions at the Jet 
Propulsion Laboratory (JPL). 

How accurately can the thermal recoil force be deter- 
mined? It always has been recognized that accurate com- 
putation of this force is a difficult task. Computing the 
total amount of heat emitted by a spacecraft is straight- 
forward: if the thermal power of internal heat sources is 
known and the spacecraft is in a steady state, the amount 
of heat it radiates must equal the generated heat. 

The recoil force, however, depends not only on the to- 
tal amount of heat radiated by the spacecraft, but on 
the differences in heat radiated in different directions. 
These, in turn, are calculated using detailed knowledge of 
the spacecraft's geometry and material properties, which 



® For instance, Thermal Desktop, the Thermal Radiation Analysis 
System (TRASYS) and Systems Improved Numerical Differenc- 
ing Analyzer (SINDA) [rlliot. 
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may be poorly known. 

Considering the force model we note that once 

a comprehensive numerical model is available, it needs 
to be evaluated only a modest number of times in order 
to obtain the efficiency factors through Eq. (|47)) . Ad- 
ditional evaluations can be used to determine the error 
of this homogeneous and linear approximation, which we 
denote with (Tmodci- If we assume that sources of error 
are independent and not correlated, the error of the recoil 
force estimate (|46l) can then be written as 

1 " 

''Fit) - - ^ + ^'^B. (t)) + aLdeb (48) 

1=1 

where uncertainties in the knowledge of Bi{t) are rep- 
resented by the standard deviations aSiit), and uncer- 
tainties in the calculation of are represented by the 
standard deviations cr^. . 

For spacecraft, the values oi Bi[t) are either measured 
values from telemetry, or nominal values from design doc- 
umentation. In the former case, ctb. {t) can be obtained 
by considering sensor sensitivity and telemetry resolu- 
tion. In the latter case, uncertainties may be available 
from design documentation, or may be inferred, e.g., by 
comparing nominal power consumption with measured 
values of available electrical power. 

The values of the cr^. are the most difficult to estimate. 
These values must be developed on a case-by-case ba- 
sis, accounting for uncertainties in the knowledge of the 
spacecraft's geometry, material properties, and physical 
configuration, as well as the effects of aging on these. 

VI. ORBIT DETERMINATION 

The position of a distant spacecraft is rarely ob- 
served directly. Instead, the spacecraft's position is 
inferred from radio-metric observables, notably radio- 
metric Doppler and range measurements. The expected 
values of these measurements can be computed if the 
spacecraft's position and velocity are known. Not in- 
cluding general relativistic corrections [13, [l3| , the space- 
craft's equation of motion can be written in the general 
form 

-E^^^.^ + ^E^- (49) 

where r is the position of the spacecraft, G is Newton's 
gravitational constant, Mj are the masses and are the 
positions of bodies that influence the spacecraft's position 
gravitationally (these would typically include all major 
solar system bodies, and any smaller bodies that are suf- 
ficiently near the spacecraft to affect its orbit), m is the 
spacecraft's mass, and are any nongravitational forces 
acting on the spacecraft. 

Once the orbit of the spacecraft is known, the expected 
values of radio-metric observables can be calculated by 



taking into account signal propagation in the solar system 
environment, atmospheric effects, antenna locations, the 
relative motion of antennae and spacecraft, and other 
effects such as the spacecraft's spin that can influence 
the signal. 

Equation is a second-order differential equation in 
r. To obtain a specific solution of such an equation, one 
needs suitable initial conditions, which can be specified 
in the form of the initial state vector that consists of 
r(to) and r(to) at some time t = to. The purpose of 
the orbit determination effort is to find initial conditions 
for which the difference between computed and observed 
radio-metric values, i.e., the residual, is minimized. 

The small forces model in (|49l) may be parameterized. 
If the parameter values are unknown, the orbit determi- 
nation exercise can be used to find these along with the 
initial state vector. 

As an example, solar pressure may be represented as a 
force acting on the spacecraft, modeling the spacecraft's 
geometry, with given solar absorptance, specular, and 
diffuse refiectivity coefficients of its surfaces. These co- 
efficients may not be known in advance accurately, but 
they can be determined along with the orbital initial con- 
ditions through minimizing the residuals of radio-metric 
observables by simultaneous adjustment of the initial 
state vector and the solar pressure coefficients. 

Precision orbit determination requires the repeated 
evaluation of ((49|) a large number of times. Therefore, 
it is essential that the forces F^ can be computed with 
as little computational overhead as possible. In the case 
of the thermal recoil force, this precludes the possibility 
of evaluating a comprehensive thermal model multiple 
times at every point along the spacecraft's orbit. How- 
ever, after the coefficients have been determined by 
evaluating the comprehensive thermal model a limited 
number of times, the thermal recoil force ([45]) can be 
incorporated easily into the orbit determination process. 

In a more innovative approach, we can treat the co- 
efficients as unknown parameters, and use the orbit 
determination program to determine their values, along 
with the initial state vector and other parameters. This 
approach is especially notable because it is applicable 
even when no detailed thermal model for a spacecraft 
is available. When a detailed thermal model is present, 
agreement between the two methods validates that model 
and the hypothesis that no additional forces of unknown 
origin act on the spacecraft. Conversely, if the two cal- 
culations are in significant disagreement, that can be a 
strong indication that an anomalous force of unknown 
origin is present. 

So long as the amount of heat generated by on-board 
components is known and the spacecraft's geometric and 
optical properties are not changing with time, one can 
make the assumption that the on-board generated recoil 
force is well modeled by (|46)) and that no other unknown 
forces affect the spacecraft, and then proceed to verify 
this model and at the same time determine the values of 

by fitting to the orbital data. Conversely, even when 
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a detailed thermal model is available, independent deter- 
mination of the from orbital data alone offers a robust 
way to verify the results of the thermal model and help 
confirm or reject any hypotheses concerning the presence 
of additional forces of unknown origin. 



VII. THE CASE OF THE PIONEER ANOMALY 

We have applied the methods that we developed in the 
previous sections to the study of the Pioneer anomaly [1] . 

Pioneer 10 and 11, launched in March 1972 and April 
1973, respectively, were the first man-made objects to 
travel to the outer regions of the solar system and be- 
yond. After flying by Jupiter (and, in the case of Pi- 
oneer 11, Saturn), the spacecraft continued on hyper- 
bolic escape trajectories, while they were being tracked 
by NASA's Deep Space Network system of radio track- 
ing stations. Pioneer 11 remained operational until 1995, 
although precision navigation of this spacecraft ended in 
1990 due to on-board failures. Pioneer 10 was operat- 
ing as late as 2003, and precision Doppler measurements 
were received from this spacecraft until the end of its 
mission'^. 

The Pioneer spacecraft [l^ were spin stabilized. Their 
spin axis, coinciding with the antenna axis, was pointed 
towards the Earth to ensure continuous communication. 
Spin stabilization meant that trajectory correction ma- 
neuvers were infrequent; most of the time late in their 
missions. Pioneer 10 and 11 were flying undisturbed. Be- 
cause of this, the twin Pioneers were considered a reli- 
able platform for precision gravitational measurements, 
searching for a planet beyond Pluto, and for gravitational 
waves. 

While neither "planet X" nor gravitational waves were 
detected, there remained an unexplained residual be- 
tween calculated and observed Doppler data 0, Q . This 
residual can be eliminated by assuming that a constant 
acceleration of unknown origin pushes the spacecraft to- 
wards the Sun. The magnitude of this acceleration is 



ap = (8.74 ± 1.33) x 10"'" m/s 



(50) 



A possible origin of this approximately constant accel- 
eration could be a thermal recoil force. Electrical power 
on board the Pioneer spacecraft was generated by low- 
efficiency RTGs that produced waste heat in excess of 
2 kW [1[ throughout the mission. Electrical equipment 
on board the spacecraft produced an additional ~ 100 W 
of heat (see Fig. [T|). This heat was radiated away by the 
spacecraft in a complex pattern, as determined by the 
spacecraft's geometry and material properties. A small 
anisotropy, less than 2% in magnitude, would be suffi- 
cient to provide the necessary force to yield the anoma- 
lous acceleration dST 
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FIG. 1: RTG heat (red -I-) and electrical heat (green x), 
measured in W, on board Pioneer 10, from telemetry [H, 



The RTGs are mounted on the Pioneer spacecraft at 
the end of booms that are approximately 3 m in length. 
The RTGs are compact objects. This suggests that the 
recoil force due to RTG heat may be modeled accurately 
as a homogeneous linear function of the RTG heat, in 
accordance with the discussion in Sec. |TVl Further, al- 
though each spacecraft has four RTGs, they are mounted 
symmetrically and their temporal behavior is nearly iden- 
tical: therefore, they can be treated as a single heat 
source. 

Electrical heat is produced inside the spacecraft body. 
Most of this body is covered by multilayer thermal insula- 
tion, resulting in small exterior temperature differences. 
This implies that the recoil force due to electrically gen- 
erated heat can also be a homogeneous linear function of 
the electrical heat, and further, the possibility that any 
particular distribution of heat sources inside the space- 
craft body can be neglected, only their total thermal out- 
put must be considered. 

Additional heat sources on board the spacecraft in- 
clude 12 small (1 W) radioisotope heater units (RHUs) 
and the propulsion system. The former can be ignored 
due to their geometry; most of the heat produced by the 
RHUs is radiated in a direction perpendicular to the spin 
axis^. Heat generated by the propulsion system, in turn, 
can be ignored as these events are transient, and any 
thermal recoil force due to propulsion system heating is 
masked by uncertainties in the modeling of the maneu- 
vers. 

Further, as the Pioneer spacecraft are spinning, we 
only need to consider the recoil force in the spin axis 
direction, according to Sec. IIVI 

Temperatures inside the spacecraft are nearly con- 
stant, changing only on the times cale of years. There- 
fore, the spacecraft is accurately described by a steady- 
state model. 

Using recently recovered documentation 1], a highly 
detailed finite element model incorporating '--^3000 nodes 
and 2600 plate elements, using 3.4 x 10^ radiation conduc- 
tors and ^7000 linear conductors has been constructed 



Only Doppler; Pioneer 10 and 11 had no range observable. 



* Private communication from Jim Moses, TRW retiree. 
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[7|. Temperatures and power readings from telemetry 
were used as a redundant data set to establish boundary 
conditions. Analysis of this model is presently under way 
and results will be reported when they become available. 

This model will also be used to validate the assump- 
tions leading to (|46|) . notably by verifying that the result- 
ing recoil force is indeed a homogeneous linear functional 
of the RTG and electrical heat: 

F(<) = -{^rBrit)+^eBe{t))s, (51) 
C 

where and are the efficiency factors associated with 
RTG thermal power Br{t) and electrical power Be{t). 
As before, s is a unit vector pointing in the spin axis 
direction. 

The discrepancy between the linear model (|5ip and a 
comprehensive thermal model can be evaluated to yield 
an estimate of crmodoi. We have not yet computed the 
Cmodoi numerically. 

The value of Be{t) is available from telemetry (Fig.[T]). 
Uncertainties in this value are due primarily to two fac- 
tors. First, telemetry has limited resolution (analog sen- 
sor readings are telemetered after conversion to 6-bit bi- 
nary values). Second, for some instruments on board, 
only their nominal power consumption is known, teleme- 
try provides only their on/off state, not their actual 
power level. Taking all these uncertainties into account, 
we calculate 

ag^ = 1.8 W. (52) 

The total power of the RTGs is known precisely from 
prelaunch documentation. The physics of the radioactive 
decay of the ^'^®Pu fuel is well understood. The amount 
of power removed from the RTGs in the form of electrical 
energy is telemetered to the ground, and thus the time 
dependence of the RTG thermal power Br{t) is known. 
The primary source of uncertainty in the calculation of 
Br{t) is the limited resolution of this telemetry. This 
uncertainty is calculated as 

o-B, = 1.1 W. (53) 

We note that as^ and ctb^ are anticorrelated; if the 
electrical heat is overestimated, the RTG heat is underes- 
timated using the same telemetry, and vice versa. Treat- 
ing the two sources of error as uncorrelated, therefore, 
yields a conservative error estimate. 

Additional temperature information is available in the 
telemetry stream, measured by sensors located at various 
points around the spacecraft. These temperature read- 
ings offer redundant information about the thermal state 
of the spacecraft that can be used to verify and validate 
thermal models. 

The effort to develop a comprehensive thermal model 
of the Pioneer 10 and 11 spacecraft is on-going. This 
effort is expected to determine an estimate for and ^r- 
A naive analytical model of the spacecraft, verified by 



a simple ray-tracing computational model, suggests that 
the values must be approximately 

=i 0.36, ^ 0.010, (54) 

albeit the relative error on these figures may be a high as 
100% or more, due to the simplicity of the models that 
were used to obtain them. Nevertheless, we use these 
figures as typical figures in the present analysis. 

We have recentl y d eveloped a precision orbit deter- 
mination program [l9| that can process Pioneer 10 and 
11 Doppler data. It uses the latest JPL ephemerides 
(DE-414) to determine the position of solar system bod- 
ies. The program models spacecraft orbits using rela- 
tivistic equations of motion. It accurately models sig- 
nal propagation by taking into account, for instance, the 
Shapiro time delay and effects of the troposphere and so- 
lar plasma on the radio signal. This program has been 
used successfully to confirm the existence of the Pioneer 
anomaly . The program also has the capability to uti- 
lize spacecraft telemetry and model on-board generated 
thermal recoil forces 

An effort to recover all available Pioneer 10 and 11 
data is presently on-going [l], 0]. Before we apply our 
method to this soon complete data set, it was essential 
to demonstrate the viability of our method. Notably, we 
would like to know if it is possible, in principle, to dis- 
tinguish between a constant sunward acceleration and a 
thermal recoil force. For this purpose, we built simulated 
sets of Pioneer 10 Doppler data. In one particular simula- 
tion, we used the following fictitious values of a constant 
acceleration term oq and thermal coefficients and S,r'. 

ao = 2x IQ-^" m/s^ (55) 
= 0.3, (56) 
= 0.015, (57) 

consistent with (jM]) . 

Furthermore, the simulation utilized actual Pioneer 10 
telemetry to model the internal heat of the spacecraft. 
The simulated data set ran from 1987 to 1998 and con- 
tained 13,534 Doppler data points. To make the simu- 
lation realistic, Gaussian random noise with ct = 5 mHz 
was added to the Doppler data. Additionally, a sinu- 
soidal diurnal term and a sinusoidal annual term, both 
with a peak-to-peak amplitude of 10 mHz, were added to 
the signal, to simulate possible mismodeling, by effects 
such as those of the atmosphere on the signal and of the 
dynamics of the solar system. 

Position values rounded to the nearest 1000 km and 
velocities rounded to the nearest m/s were used as the 
initial state vector (this is the typical magnitude of er- 
ror we observe when we use Pioneer ephemeris data from 
JPL Horizons On-Line Ephemeris System^ for initial con- 
ditions.) We used initial values of oq = 1 x 10^^" m/s^. 



^ |http : //ssd ■ jpl ■ nasa . gov/?horlzons | 
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FIG. 2: Simulated Pioneer 10 Doppler prefit residuals. 
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FIG. 3: Post-fit residuals of the simulated Pioneer 10 data set 
after the values of ao, (,e, and along with the initial state 
vector were fitted successfully. 

= 0, and = 0. Though we have the capability to 
deal with maneuvers in the actual data, for the purposes 
of this exercise we did not simulate maneuvers. The re- 
sulting prefit residuals are shown in Fig. [2] 

The goal of this simulation was to demonstrate that 
even in the presence of noisy data, a constant acceleration 
term and acceleration due to thermal radiation can be 
clearly distinguished. Despite the presence of noise, our 
orbit determination algorithm successfully recovered the 
values of 

ao = (2.1107 ± 0.0170) X 10^1° m/s^ (58) 
= 0.292 50 ± 0.002 54, (59) 
= 0.014 856 ±0.000 081. (60) 

The post-fit residuals are shown in Fig. [31 The root mean 
square residual of this solution is 5.84 mHz, which cor- 
responds to the noise that was added to this simulated 
data set. 

Furthermore, the cross-correlation between oq, ^e, and 



TABLE L Covariance matrix elements for ao, and ^e- 





ao 




& 




ao 


2.91 X 10" 


24 


-2.15 X 10"" 


-2.95 X lO-^'^ 


& 


-2.15 X 10' 


-19 


6.51 X 10"^ 


-1.20 X 10"^ 




-2.95 X 10" 


-15 


-1.20 X lO"'^ 


6.45 X 10"^ 



remains small. This can be seen by visual inspection of 
the relevant elements of the covariance matrix, shown in 
Table m We note that, after normalizing using the values 
of flo, S,r, and ^e, the diagonal elements of the covariance 
matrix dominate. 

These results indicate that the approach we presented 
is feasible. It is possible, in principle, to distinguish a 
constant acceleration from time- varying acceleration due 
to thermal radiation using Doppler data alone, even when 
the data has a moderate amount of noise. Nonetheless, it 
is imperative to reduce the noise in the data as much as 
possible, for example by carefully modeling small effects 
such as those of the atmosphere and solar plasma on the 
spacecraft's radio signal, or small accelerations due to 
fuel leaks and maneuver uncertainties. 

Once an improved thermal model becomes available, it 
can be used to verify the linear hypothesis expressed in 
(I46p . which forms the basis of the approach we present 
here. The thermal model may also be used to quantify 
the error margins on and • On the other hand, anal- 
ysis of recently recovered Doppler data can confirm if the 
orbital behavior of the Pioneer spacecraft remained con- 
sistent throughout their missions, and may also help re- 
duce the error margins on any residual acceleration that 
remains after accounting for the thermal recoil force. 

VIII. CONCLUSIONS 

An object that emits heat experiences a recoil force 
due to radiation pressure. In this paper, we developed 
the basic equations that can be used to estimate the mag- 
nitude of this recoil force, and relate the recoil force to 
the amount of heat produced internally. We have been 
able to show how, under specific circumstances, the recoil 
force can be modeled as an homogeneous linear function 
of the power of discrete internal power sources. When 
this approach is applicable, the linear relationship can 
be readily incorporated into orbit determination efforts. 

To analyze the trajectory of Pioneer 10 and 11, we de- 
veloped orbit determination software that estimates the 
thermal recoil force acting on the spacecraft. Our soft- 
ware uses telemetry information as it calculates the ther- 
mal power of on-board heat sources as functions of time. 

A comprehensive thermal model, presently under de- 
velopment, will allow us to verify the key assumptions 
behind our modeling, most notably the assertion that the 
thermal recoil force is accurately modeled as a linear, ho- 
mogeneous function of electrical heat and heat from the 
radioisotope thermoelectric generators. 
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Using a simulated Doppler data set and actual Pioneer 
10 telemetry, we demonstrated that it is possible in prin- 
ciple to distinguish acceleration due thermal radiation 
from a constant sunward acceleration term. 

Newly recovered Doppler data are now available as a 
result of an extensive data recovery effort [J, 0] ■ This will 
allow us to extend our analysis, and verify whether or not 
the thermal recoil force can account for the anomalous 
acceleration of Pioneer 10 and 11. These results will be 
published elsewhere when they become available. 

We emphasize that the approach presented here, no- 
tably the direct utilization of flight telemetry in precision 
spacecraft navigation codes, has never been attempted 
before. The approach we describe is applicable not only 
to the case of Pioneer 10 and 11, but also to the case of 
present and future spacecraft. One mission in particular 
that may benefit from this approach is New Horizons, on 
its way towards an encounter with Pluto in 2015. While 
presently not used for gravitational research, such inves- 
tigations could be conducted during its multiyear cruise. 
If such an investigation is undertaken, it will require ac- 
curate estimates of the thermal recoil force due to the 
waste heat produced by New Horizons' RTG and electri- 



cal equipment. 
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